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Abstract. We study the halo bispectrum from non-Gaussian initial conditions. 
Based on a set of large A-body simulations starting from initial density fields with local 
type non-Gaussianity, we find that the halo bispectrum exhibits a strong dependence 
on the shape and scale of Fourier space triangles near squeezed configurations at large 
scales. The amplitude of the halo bispectrum roughly scales as /^ L . The resultant 
scaling on the triangular shape is consistent with that predicted by Jeong & Komatsu 
based on perturbation theory. We systematically investigate this dependence with 
varying redshifts and halo mass thresholds. It is shown that the /nl dependence of 
the halo bispectrum is stronger for more massive haloes at higher redshifts. This 
feature can be a useful discriminator of inflation scenarios in future deep and wide 
galaxy redshift surveys. 
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1. Introduction 

The standard cosmological model has successfully explained the observed statistical 
properties of the cosmic microwave background (CMB) radiation and the large scale 
structure (LSS) traced by galaxies (e.g., [1, 2]). The model usually assumes that 
the primordial density/temperature/curvature fluctuations follow Gaussian statistics. 
Recently, however, possible deviations from standard Gaussian statistics has attracted 
great attention with rapid progress in observational techniques. It offers an opportunity 
to access cosmological information beyond traditional power spectrum analysis. Many 
recent works have discussed the constraints and future detectability of possible 
deviations from Gaussianity through the observations of CMB and LSS (e.g., [3, 4]). 

According to the inflationary scenarios, primordial curvature perturbations are 
generated during the accelerated phase of cosmic expansion. The simplest inflation 
models, in which the inflation takes place with the slow-roll single scalar field that has 
a canonical kinetic structure, generally predicts a nearly scale-invariant spectrum of 
curvature perturbations, and a small departure from Gaussianity. On the other hand, 
a variety of inflation models that produce a large non-Gaussianity has been recently 
proposed (see [5] for a review). Among these, the models with large non-Gaussianity 
generated by non-linear dynamics on super-horizon scales can have a generic prediction 
for non-Gaussian properties. Denoting a Gaussian field by <&q(x), the Bardeen's 
curvature perturbation during the matter era is characterized as [3]: 



This type of non-Gaussianity, described as a local function of the Gaussian field, is often 
called local type non-Gaussianity, and the leading-order coefficient /nl, which controls 
the strength of non-Gaussianity, has information on the generation mechanisms for non- 
Gaussian fluctuations. Although the current constraint on the parameter /nl from CMB 
data is —9 < /nl < 111 (95%C.L.) [1] and it is still consistent with Gaussian (/nl = 0), 
the constraint will be tightened by the on-going CMB experiment, Planck [6]. As 
standard inflation models generally predict |/ NL | <C 1, detection of large non-Gaussianity 
immediately implies the non-standard mechanism for generation of primordial curvature 
perturbations. 

In this paper, we focus on how this type of non-Gaussianity alters the statistical 
properties of LSS. Let us first define the power spectrum of mass density fluctuations 
assuming isotropy and homogeneity: 



where S m (k) denotes the Fourier transform of the density contrast, while 5j^(k) 
represents the Dirac delta function. If the density field follows Gaussian statistics, its 
power spectrum determines all the statistical properties. Next we define the bispectrum 
of a mass density field: 



$(*) = * G (x) + /nl {& g {x) - (& G (x))} + ■■■. 



(1) 




(2) 



(5 m (ki)5 m (k 2 )5 m (k 3 )) = (2tt) 



3 <5 D (fci + k 2 + fe 3 )-B m (fei, fc 2 , fe 3 )- 



(3) 
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Since this is the lowest-order non- vanishing quantity in the presence of non-Gaussianity, 
the bispectrum is naively expected as the best measure for non-Gaussianity (e.g., 
[7, 8, 9, 10]). 

Recently, however, the galaxy or halo power spectrum has been reconsidered in the 
presence of local and/or equilateral type primordial non-Gaussianity both analytically 
and numerically (e.g., [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]). The matter 
power spectrum or bispectrum is not a direct observable, and the real measurement 
of LSS gives galaxy power spectrum or bispectrum, as defined similarly to equations 
(2) and (3). Since galaxies are biased tracers of the dark matter distribution, the 
information of /nl is imprinted in a different manner: a new contribution coming from 
the primordial non-Gaussianity may dominate over the Gaussian term in the galaxy 
power spectrum, P g {k), at very large scales (k ^ 0.01/iMpc -1 ). This new contribution, 
which is sometimes referred to as the "scale-dependent bias" of the galaxy power 
spectrum, may be a powerful indicator to constrain / NL . Indeed, it has been recently 
applied to the clustering statistics of SDSS LRG and quasar samples, and the tight 
constraints on /nl are comparable to those obtained from CMB measurements have 
been obtained [19]. 

The purpose of this paper is to examine the bispectrum of biased tracers in detail. 
While the matter bispectrum in the presence of primordial non-Gaussianity has been 
studied in the literature using both perturbation theory and numerical simulations, the 
galaxy bispectrum may significantly differ from the matter bispectrum in the presence of 
primordial non-Gaussianity, just like the difference in the power spectra. Since the local 
type non-Gaussianity can be straightforwardly implemented within iV-body simulations, 
numerical study on the bispectrum for the dark matter haloes is the first important step 
toward a practical understanding of the galaxy bispectrum. 

Incidentally, Jeong and Komatsu (2009) recently proposed a new parametrized 
model for the halo/galaxy bispectrum ([23], see also [21]) based on the peak bias model 
[24] and the local bias model [25] . They found that the formula for the galaxy bispectrum 
used in [10] was missing important contributions from the scale dependent bias effects 
and they discovered new terms that are important at "squeezed" configurations where 
fci,/c2 ^> k 3 . It was argued that these new contributions enable us to put stronger 
constraints on /nl than those obtained in [10]. It is of great importance to confirm the 
scale-dependent bias effects in the bispectrum by iV-body simulations. 

The rest of this paper is organized as follows: we first review the analytical models of 
the power spectrum and the bispectrum of biased tracers in section 2. We then describe 
the setup and initial conditions for iV-body simulations in section 3. As a first check 
of our simulations, in section 4, we compute the matter and halo power spectra, and 
the results are compared with previous works. Section 5 gives the main results of this 
paper, in which the simulation results for the matter and halo bispectra are presented 
and compared with predictions from analytic models, particularly focusing on their 
scale dependence. The dependence of the halo bispectrum on the halo mass threshold 
and redshift is also investigated in detail. Section 6 discusses the future prospects for 
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detecting the primordial non-Gaussianity using the scale-dependent properties of the 
halo/galaxy bispectrum. Finally, section 7 is devoted to conclusions and discussions. 

2. Theoretical models 

In this section, we summarize the theoretical predictions of the power spectrum and 
bispectrum. We use perturbation theory to examine the matter power spectrum and 
bispectrum, and then present those of biased tracers based on the local bias model. 
For the scales of our interest (k £ O.lh Mpc" 1 ), the non-linearity of gravitational 
evolution is moderate and the perturbation theory is valid and trustful. We especially 
focus on the behavior of the bispectrum in the squeezed limit on large scales, where 
k = ki = k 2 = ak 3 , k — » and a>l. 

Let us first consider the matter density fluctuation. We perturbatively expand this 

as 

5 m (k; z) = 5«(fc; z) + 5%\k; z) + 5^(k; z) + ■ - ■ . (4) 

The linear-order solution is related to the Bardeen's curvature perturbation in equation 
(1) in Fourier space by 

6^(k;z)=M(k;zMk) } (5) 

where the conversion factor A4 is defined as 

In the above, Q m is the current matter density normalized by the critical density, H is 
current Hubble constant, T(k) denotes the matter transfer function normalized to unity 
at k — > and D(z) is the linear growth rate normalized to the scale factor in the limit 
of matter dominant era. The higher-order solutions, 8j£'(k] z), are formally written as 

S^(k; z) = J d 3 Ql ■ ■■d 3 q n 6 D (k - qi ... n )F n ( qi , . . . , q n )S^(q i; z) ■ z), (7) 

where qr x ra = q x + • • • + q n , and F n are the kernel functions (see [26] for a review). 

Then keeping terms up to fourth order in 5$, the power spectrum and bispectrum of 
the matter density fluctuations are given by (e.g., [13, 10]) 

P m (k; z) = Po(k; z) + 2 J tS^ 2 (<7, k - q)B (-k, q,k - q, z) 



+ 2 J 0- i {F 2 (q 1 k- q )} 2 P o (q;z)P o (\k-q\;z) 

/d?pd? q 
(2?r) 6 F2 ^ p ' k ~ p^ 2 ^' ~ g ) T °(P' k ~ P> ~ k ~ 9; z), 

+ Po(k] z) J S^Fs(k, q, -q)P (q; z), (8) 



B m (ki, k 2 , h; z) = 2f] 



NL 



P {k 1 ;z)P {k 2 ;z)M{k 3 ;z) , 

- + (cyc y 



M(h; z)M(k 2 ; z) 
+ 2F 2 (k 1 , k^Poih; z)P (k 2 ; z) + (eye), (9) 
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where (eye.) denotes the cyclic permutations over the indices and P , B , and T are 
the power-, bi-, and tri-spectra of 5$. In equation (8), the second term is the first non- 
trivial correction in the presence of primordial non-Gaussianity, and the function Bo 
implies the primordial bispectrum, which corresponds to the leading-order contribution 
in equation (9). Note that the contribution coming from the primordial trispectrum To 
is small for local type non-Gaussianity with reasonable values of /nl, and we drop this 
term in computing P m (see [13]). 

On the other hand, the power spectrum of biased tracers has been recently discussed 
in the literature [11, 14, 19, 12, 20, 13], based on several different formalisms including 
peak bias, halo bias according to the peak-background split, and local bias. The 
resultant expressions of the galaxy/halo power spectrum are basically the same, and 
are summarized in the form 

P g (k; z) = b\ jl + 2 | / NL M-\k- z) J P (k; z), (10) 

where b\ and 62 are the bias parameters relating the galaxy overdensity to the matter 
overdensity. In Appendix A, we present a derivation of (10) based on the local bias 
formalism. The explicit expressions for the bias parameters h\ and 62 can be obtained 
both from the peak bias and halo bias formalisms, and their results are basically the 
same in the high-peak/threshold limit. In Appendix B, we show that in the high-peak 
limit, there is a clear relationship between the peak bias and the local bias prescriptions, 
and the parameters bi and b% are related to each other in terms of the critical density 
of the spherical collapse model, 5 C ~ 1.686, as 6 2 = <5 C (61 — 1). 

In equation (10), the factor in the braces manifestly depends on the scale, which 
is the main source for "scale-dependent bias". On large scales (k — > 0), the function 
M. is roughly proportional to k 2 , and it strongly affects the galaxy power spectrum. 
In this respect, the scale-dependent property will be a clear indicator of primordial 
non-Gaussianity of the local type, and it has been extensively tested against N- 
body simulations. Several recent studies have suggested that some modifications to 
this formula are required in order to model the scale-dependent bias more accurately 
[15, 16, 18]. For example, [18] proposed a slight modification to the relation of bias 
parameters, which reproduces results from iV-body simulations very well: 

62 = ^5(61-!), (11) 
with q = 0.75, which comes from the ellipsoidal collapse model. 

Now, we turn our focus to the bispectrum of biased tracers in the presence of 
primordial non-Gaussianity. According to the analytical study by [23] (see also [21]), the 
effect of local-type primordial non-Gaussianity is mainly imprinted on the bispectrum 
of squeezed triangular configurations. The galaxy bispectrum is then expressed as 

B g {k, a; z) = Bf» (k, a; z) + / NL (k, a; z) + f 2 L B® (k, a; z), (12) 

where each term of the right-hand side of this equation has the following asymptotic 
form: 

Bf\k, a- z) ~ bj 62 PL {z) k 2n ° a , (13) 
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( y + l(k,a;R) 



2», 



" 2 a\ 



8^6aitM)(«)^i^o(«)* 



2(i, s 



" 4 a 3 . 



(14) 
(15) 



Here, we focused on the isosceles triangles, and parametrized their dependence as 
k = ki — k 2 = a/c 3 , and n s denotes the scalar spectral index. See Appendix A for 
a more rigorous expression. The function I(k, a; R) weakly depends on k, a and the 
smoothing scale, R, and it can be approximated as 26/7 + X(k,a; R) ~ 34 on large 
scales. In the above, we take the limit of k — > and a ^> 1, and -Pfc^O) -Mfc->-o are 
defined through P (k; z) — > P k ^ (k)k ns , M(k; z) — > Mk^o{z)k 2 . See Appendix A and 
also [23] for more details. 

Similar to the scale dependence of the galaxy/halo power spectrum in equation 
(10), the amplitude of bispectrum is also affected by the primordial non-Gaussianity 
in a scale-dependent way. In particular, the term which is of quadratic order in 
/nL) becomes dominant on large scales and exhibits a behavior strongly dependent on 
k and a. Thus, we might conclude that this is the most important indicator of /nl, 
capable of detecting the primordial non-Gaussianity with upcoming galaxy surveys of 
large volumes. 

In what follows, we will examine this scale dependence in the bispectrum of 
simulated dark matter haloes, with particular attention on the squeezed configurations 
on large scales. 



3. iV-body Simulations 

3.1. Setup 

We adopt the WMAP5 best-fit flat ACDM model [1]: O m = 0.279, tt A = 0.721, 
f2 b = 0.046, h = 0.701, a 8 = 0.817 and n s = 0.96, where Q m , Q\ and 1\ are the 
matter density, cosmological constant and the baryon density normalized by the critical 
density, h is Hubble constant normalized by 100kms~ 1 Mpc~ 1 , ag is the r.m.s. linear 
density fluctuation smoothed by a top hat window function with radius of 8/i~ 1 Mpc and 
n s is the scalar spectral index. We calculate the linear matter transfer function using the 
CAMB code [27]. We have completed a total of 140 realizations of matter clustering data: 
20 per each of the seven values for the local-type primordial non-Gaussianity parameter, 
/nl = 0, ±100, ±300 and ±1000. All the simulations were ran with the Gadget2 code 
[28]. We adopt N = 512 3 particles in boxes of side 2000/i _1 Mpc, and set the softening 
length being 0.2/i _1 Mpc. The parameters adopted in the simulations are the same as in 
[29, 30]. We have tested the mass/force resolution by changing the box size and found 
that the results are well converged at large scales (~ 1% accuracy at k £ 0.3/iMpc" 1 ) 
for the matter power spectrum. 

In setting the initial conditions, we first generate a random Gaussian field, $Gj 
whose power spectrum is proportional to k nB ~ A . We then apply the inverse Fourier 
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transform, and add the non-Gaussian contributions in real space according to equation 
(1). Finally the real-space quantity is transformed back to the Fourier space, and 
converted to the linear density fluctuations by multiplying J\4(k; z) defined in equation 
(6). We use second-order Lagrangian perturbation theory to calculate the displacement 
field for 512 3 particles placed on a regular lattice [31] at z = 31. 

We store outputs at z — 2, 1 and 0.5, and identify haloes for each output using 
a FOF group finder with linking length of 0.2 times the mean separation. We select 
haloes in which the number of particles, N, is equal to or larger than 10, corresponding 
to the haloes with masses 4.6 x 1O 13 /i _1 M . We also analyze haloes with N > 20 and 
N > 30 to see the dependence on halo mass. 

As a first check of the reliability of our simulations, we show in Fig. 1 the ratio 
of halo mass function with and without primordial non-Gaussianity, R(M; /nl) = 
n n G(M; fm,)/nG(M), at z = 0.5, where different colors correspond to different values of 
/nl: 1000, 300, 100, 0, -100, -300 and -1000 from top to bottom. Both simulations 
and theoretical models suggest that the local type non-Gaussianity alters the mass 
function at the high-mass tail in the literature: a positive (negative) /nl enhances 
(suppresses) the tail. As shown in Fig. 1, our simulations agree well with previously 
proposed analytical models. The plotted lines show the model proposed by [32] (solid): 



R(M; / NL ) = 1 + 

o o. 



cc 




CC 



(16) 



based on the Edgeworth expansion to the probability density function, and by [33] 
(dashed): 



R{M; /nl) = exp 





s 3 \ 




0-i) 



1 5 Pr dS* 



6Jl-5 cc S 3 /3dlna R 



\Jl- SccSz/S 



;i7) 



obtained by the saddle-point approximation to the level excursion probability. In the 
above, ctr = (S 2 ^) 1 ^ 2 and S 3 = (Sf-^/a^. These quantities are given as the function of 
mass M through the relation M = (47r/3)p m R 3 , and linearly extrapolated to z = 0.5. 
We also define 5 ec = qS c again motivated by ellipsoidal collapse model [18] . See also [16] 
for another model designed to fit to their iV-body simulations. 

There are several claims on the systematics for the estimated halo mass by FOF, 
and in fact our mass function fits better with [33] when we correct that effect using the 
empirical formula of [34]. We conclude here that our halo catalog is accurate enough to 
investigate its clustering statistics. 



3.2. Measurements of the power spectrum and bispectrum 

Here, we briefly mention how to measure the power spectrum and bispectrum in our 
simulations. 

We assign particles (or haloes) to 1024 3 grid points using the Cloud- in- Cells (CIC) 
algorithm [35]. We then Fourier transform the density field, and divide each mode by 
the Fourier transform of the CIC kernel to correct for the effect of assignment. We 
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M[h- 1 M sun ] 

Figure 1. The ratio of the halo mass functions for non-Gaussian and Gaussian initial 
conditions at z — 0.5. The symbols show the measurements from our simulations, 
while the lines are equations (16) and (17). The values of /nl are 1000, 300, 100, 0, 
-100, -300, and -1000 from top to bottom. 

made sure that the results are well converged at the scale of our interest by changing 
the number of grid points. We logarithmically divide the measured power spectrum 
and bispectrum into wave number bins starting from /c min = 0.003/iMpc" 1 and with 10 
bins per decade. We select "isosceles" triangles, whose two longer sides, ki and &2, fall 
into the same bin for the bispectrum analysis. In this sense our isosceles triangles are 
not strictly isosceles, but we adopt this convention to reduce the statistical errors on 
the measured bispectrum caused by small number of triangles in /c-space. In plotting 
results, we assign each data point to the logarithmic-central value of wave number in 
that bin. 

4. Results of Power Spectrum 

In this section, we present the power spectrum measured from iV-body simulations 
We first compare the measured matter power spectrum with the perturbation theory 
prediction. We then present the halo power spectrum, and compare it with the analytic 
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models proposed in the literature. These are important sanity checks to justify the 
results of our simulations in subsequent sections. 

4-1- matter power spectrum 

We first examine the matter power spectrum, in cases with non-zero /nl- Fig- 2 shows 
the fractional difference of the matter power spectra between Gaussian and non-Gaussian 
initial conditions measured at z — 0.5. The symbols represent /nl = 300, 100, 0, —100 
and —300 from top to bottom at k ^ 0.02/iMpc^ 1 . Since we use the same set of random 
seeds for the seven /nl parameters, the cosmic variance is effectively canceled out by 
taking the ratios, P^ G (/c; z)/P^(k; z). Overall, the deviations from Gaussian results 
themselves are very small (less than 1% when |/nl| = 100 in the plotted range). 

In Fig. 2, we also plot the predictions based on perturbation theory of [13], depicted 
as continuous lines. Note that recently, Ref. [36] developed another analytical model 
based on the Time-RG approach, which would be more accurate in the weakly nonlinear 
regime. However, the standard perturbation theory prediction of [13] is accurate enough 
at the scale of our interest (i.e., k £ 0.05/iMpc -1 ), as was shown by comparisons with 
TV-body simulations in Refs. [15, 16]. 

Fig. 2 shows that the results of our iV-body simulations are in reasonably good 
agreement with the model of [13], except for the case of a large non-Gaussianity with 
/nl = 300. The discrepancy between iV-body and analytic results seen in the /nl = 300 
case might be partially ascribed to the term coming from the primordial trispectrum in 
equation (8), which is neglected in the perturbation theory calculations. Nevertheless, 
the discrepancy remains at the sub-percent level, and thus does not seriously affect the 
later analysis of the matter/halo bispectrum. 

4-2. halo power spectrum 

We next consider the halo power spectrum for various values of /nl, shown in Fig. 
3. We set the minimum mass of haloes to 4.6 x 10 13 /i _1 M Q , which corresponds to 10 
iV-body particles. We also show the analytical prediction of equation (10) with the 
bias parameter bi fitted to reproduce the results of Gaussian simulations, adopting the 
relation between bias parameters in equation (11). We plot the model with q = 1 and 
q = 0.75 by dotted and solid lines, which corresponds to the original peak bias prediction 
and the fit by [18], respectively. 

Overall, the scale dependence of the halo power spectrum discussed in the literature 
can be clearly seen in our simulations with very small statistical errors, owing to the large 
total volumes. The agreement between iV-body simulations and the analytic models 
becomes better when we choose q = 0.75, consistent with [18]. 

Note, however, that the choice of q — 0.75 does not necessarily imply the best- 
fit results: q = 0.85 gives a better fit to this particular case, and the best-fit value 
of q changes with redshift and minimum halo mass. This might indicate that there 
exists some systematic effects on the halo clustering properties in our simulations. One 
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Figure 2. Fractional differences of the matter power spectra starting from non- 
Gaussian and Gaussian initial conditions. Symbols show the results of TV-body 
simulations, while lines are perturbation theory predictions of equation (8) (/nl = 
300, 100, 0, -100, -300 from top to bottom at k > 0.02/iMpc- 1 ). 



possibility is the difference in the halo finding algorithms: while we adopt the FOF 
finder, [18] use a SO finder. See also [15], where they use a FOF finder and proposed a 
fit corresponding to q = 0.8. 

We may further improve the agreement between iV-body simulations and theoretical 
predictions by including some corrections to the theory. Ref. [16] showed that the 
inclusion of two corrections coming from the changes in halo mass function and the 
matter power spectrum actually improves their results. These systematics will definitely 
be important for the application to the upcoming surveys, however, we do not pursue 
this issue in the present paper, since our primary focus is on the halo bispectrum. 

5. Results of Bispectrum 

In this section, we present the bispectrum measured from iV-body simulations. 
Throughout the analysis, we consider the isosceles triangles for the configuration of 
bispectrum, which are characterized by the two parameters k and a, defined by 
k = k\ = k% = ak 3 . We pay special attention to the squeezed triangles, a ^> 1. 
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Figure 3. The power spectrum of haloes more massive than 4.6 x 10 13 /i _1 M© at 
z = 0.5. Symbols are results of TV-body simulations, while lines are theoretical 
prediction of equation (10) with (11) where we adopt q = 1.0 for dotted lines and 
q = 0.75 for solid lines (/ NL = 300, 100, 0, -100, -300 from top to bottom). 



We first present the results of the matter bispectrum (Sec. 5.1), and then discuss how 
the halo bispectrum differs from the matter bispectrum (Sec. 5.2). While we mainly 
analyze the default halo catalog with minimum mass M min = 4.6 x 10 13 h~ l M & and 
output redshift z = 0.5, we briefly discuss how the results are changed when we vary 
the minimum halo mass and redshift (Sec. 5.2.3). 

5.1. matter bispectrum 

Let us present the results of the measured matter bispectrum. In Fig. 4, the symbols in 
each panel show the amplitude of the bispectrum measured from simulations for various 
/nl at a fixed triangle specified by k and a indicated in the panel. We also show the 
perturbation theory prediction of equation (9) by solid lines. Note that the value of a 
increases from right to left, while k increases from top to bottom. 

Although we have very large total volume, the statistical uncertainty due to 
finiteness of the simulated volume still affects the measurements. We thus take account 
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of this effect in the perturbation theory predictions: we compute the matter bispectrum 
using the second-order perturbation theory starting from linear density field realized in 
finite- volume boxes which were used to generate the initial conditions of the simulations, 
and take average over realizations. Namely, we compute 

+ ${ki)$\k a )S£\k a ) +6£\k 1 )6g\k 2 )8U(k 3 ) + 4 1) (fci)4 1) (fc 2 )4 2) (fc 3 )], (18) 

and take the average over the realizations and triangles in the bin for the perturbation 
theory prediction. As a result, the analytical predictions are in good agreements with 
measurements from simulations. 

The matter bispectrum from both simulations and perturbation theory clearly 
exhibits a linear dependence on /nl for all the configurations plotted in Fig. 4. Based 
on this results, we will discuss how the dependence on /nl is modified for the halo 
bispectrum. 



5.2. halo bispectrum 

We are now in a position to show the halo bispectrum. Since this is the first numerical 
study on the halo bispectrum in the presence of local type non-Gaussianity, it is 
important to understand the iV-body results in a model-independent manner. In this 
subsection, we first show the /nl dependence of the halo bispectrum. We then consider 
the scale dependence and compare the simulation results with the theoretical prediction 
of [23]. The dependence on the minimum halo mass and redshift is also investigated in 
section 5.2.3. 



5.2.1. /nl dependence In order to quantitatively study the /nl dependence of the halo 
bispectrum, we use all the halo catalogs with various values of /nl, and fit the measured 
bispectrum to the polynomial form: 

B h (k,a;z\h L ) =£/AiX l) (^a^), (19) 

j=0 

using the standard x 2 fitting with the variance of the data points measured from N- 
body simulations. For specific configurations of (k,a), we determine the parameters 
B^' using the halo catalogs with different values of /nl- We confirmed that the results 
are almost unchanged when we add higher-order polynomials with i > 5. 

In Fig. 5, we show the measured bispectrum as function of /nl for the same set of 
triangular configurations as plotted in Fig. 4. Note again that the value of a increases 
from right to left panels, while k increases from top to bottom panels. The fitted results 
of Eq. (19) truncating at the first order (i = 0, 1), second order (i = ~ 2) and fourth 
order {i = ~ 4) are shown respectively as dotted, dashed and solid lines. Although 
we do not show the points at /nl = ±1000 in order to focus on more realistic values of 
/nl, we take account of these results when we fit the iV-body data to Eq. (19). 
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Figure 4. The matter bispectrum. Each panel shows the results for an isosceles 
configuration specified by a = k\/k^ and k = k\ = k-2- Symbols are measurements from 
A^-body simulations (the average and the standard error among different realizations) 
and solid lines are the perturbation theory predictions of equation (9). 



The second-order term (B^ ) becomes more significant in moving from bottom 
to top and from right to left, and in the end, the top-left panel (k = 0.042/iMpc _1 , 
a = 12.6) shows strong evidence of Higher order terms {B^ and B^) seem to 

have almost no effect on the total bispectrum when |/nl| ~ 100, although they may play 
some roles at /nl = ±300 (and also ±1000). This second order term, B^ , is not seen 
in the matter bispectrum (Fig. 4) and we for the first time confirme that this really 
exists in the halo clustering in iV-body simulations. 
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Figure 5. The halo bispectrum for sonic triangular configurations. Each panel shows 
the result for an isosceles configuration specified by a = k\/kz and k = fci = &2- 
Error bars arc measurements from our simulations (the average and the standard error 
among different realizations) and solid lines are their 4-th order polynomial fits, while 
we keep the terms up to second and linear order for dashed and dotted lines. We use 
the outputs at z = 0.5 and consider the haloes more massive than 4.6 x 10 13 /i _1 Mq. 
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5.2.2. shape and scale dependence We next investigate the scale dependence of 
B^\k,a; z). We show them in Fig. 6 for /nl = 100. The left panel shows the a 
dependence when the wave number k is fixed to 0.042/iMpc -1 , while the right panel 
displays the k dependence when a = 12.6. The analytic prediction based on local bias 
(see §2 and Appendix A) predicts B^ oc k 2 a°, oc k^a 1 and B^ oc k~ 2 o? in the 
squeezed limit at large scales (k — > 0, a ^> 1), and we show these asymptotic scalings 
by short straight lines (normalizations are arbitrary). 

The a dependence in the left panel is quite consistent with the theoretical 
predictions, and the results strongly indicate that the theoretical model captures the 
nature of the shape dependence. On the other hand, the k dependence measured from 
iV-body simulations seems different from that predicted by the theoretical model. This 
implies that the wave numbers shown in the figure are not sufficiently small, and the 
approximation used in deriving the theoretical predictions is not valid. We expect 
that the asymptotic scaling appears only at the scales larger than the turn over of the 
power spectrum (i.e., T(k) ~ 1). In fact, the value of the matter transfer function at 
k = 0.042/iMpc -1 , corresponding to the wavenumber at the left-most bin in the panel, 
is T{k) = 0.293, and thus the approximation of M.{k) ~ M.k^o k 2 used in equation (12) 
cannot be applied. One can find a similar feature in Fig. 7 of Ref. [21] where the authors 
compute the galaxy bispectrum using one-loop perturbation theory adopting the local 
bias model. Although it illustrates the galaxy bispectrum for equilateral triangles, the 
asymptotic power law feature appears only at very large scale (i.e., k £ 0.03/iMpc^ 1 ). 

Nevertheless, both simulations and theory suggest that B^ > B^ > B^ at the 

(2) 

limit of small k. The B^ term will play an important role in constraining /nl from 
future surveys where we can investigate such large scales. 

5.2.3. dependence on halo mass and redshift So far, we have concentrated on the 
haloes with Mh a i > 4.6 x 10 13 /i _1 M Q at z — 0.5. In this subsection, we extend our 
analysis to the haloes with higher mass thresholds and at different redshifts to see the 
dependence of the halo bispectrum on these quantities. For this purpose, we specifically 
consider the squeezed triangle with k = 0.042/iMpc -1 and a = 12.6, corresponding to 
the configuration shown in the top-left panel of Fig. 5, and plot in Fig. 7 the amplitude 
of the bispectrum against /nl for different mass thresholds (left) and redshifts (right). 

In the left hand panel, each symbol and a line respectively correspond to the 
measurements and the polynomial fit based on equation (19) for a fixed minimum 
mass of haloes given by M min = 4.6 x 1O 13 /i _1 M (square/solid), 9.2 x 10 13 /i _1 M Q 
(triangle/dashed) and 1.4 x 10 14 h~ l M & (circle/dotted). Similar to Fig. 5, we can see a 

(2) 

clear quadratic dependence on / NL , but the role of the quadratic term B^ seems more 
significant for haloes with larger minimum masses. 

In the right hand panel, each symbol and line respectively show the measurements 
and a fit at the specific redshifts z = 0.5 (square/solid), 1 (triangle/dashed) and 2 
(circle/dotted). Here, we fix the minimum halo mass to M m i n = 4.6 x 10 13 /i _1 M Q . 
Again, one can see the quadratic dependence on / NL , which become more important for 
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Figure 6. Shape and scale dependence of the halo bispectrum. left: a dependence 
when k is fixed, right: k dependence when a is fixed. We plot the results at z = 0.5 
for haloes more massive than 4.6 x 10 13 h~ 1 MQ when /ml = 100. Error bars are 
measurements from our simulations, while dotted, dashed and solid lines arc terms 
which scale as /nlj /nL' an< ^ /nL' respectively. Short straight lines are corresponding 
analytical prediction in the squeezed limit. 



higher redshifts. 




flML fNL 

Figure 7. left: Mass dependence of the halo bispectrum at z = 0.5. Symbols and 
lines are similar to Fig. 5, but we fix the triangular configuration to be the same as 
the top-left panel, while changing the minimum halo mass: 4.6 x 10 13 , 9.2 x 10 13 and 
1.4 x 10 14 /i _1 M Q for squares, triangles, and circles, respectively. Note that although we 
do not show the results for /nl = ±1000, we take them into account for the polynomial 
fitting. Note also that the vertical axis is logarithmic above the dot-dashed horizontal 
line, while it is linear below it. right: same as left panel but the redshift dependence 
for a fixed minimum halo mass (4.6 x lQ 13 !^ 1 Mq). Squares, triangles, and circles 
correspond to z = 0.5, 1, and 2, respectively. 



Although we do not try to find the best halo catalog or optimal weighting scheme to 
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Table 1. The survey parameters used for the forecasts. We assume the one-to-one 
correspondence between haloes and galaxies to compute the minimum halo mass and 
the bias parameter from the number density. 



survey 


z 


V [/i" 3 Gpc 3 ] 


n g [/i 3 Mpc~ 3 ] 


M min [h^Mol 




IDEAL 


1.0 


100 


1 x 10~ 3 


2.8 x 10 12 


1.9 


REALISTIC 


1.0 


10 


5 x 10~ 4 


5.0 x 10 12 


2.2 


DEEP 


2.0 


3 


3 x 10" 4 


3.4 x 10 12 


3.3 



detect the signal of /nl here, the balance between having denser samplings and getting 
higher signals for /nl by selecting massive haloes is clearly very important. We will 
investigate these issues elsewhere. 

6. Prospects for future survey 

In this section, we discuss future prospects to detect /nl through the measurements of 
the bispectrum. We especially pay attention to the importance of the higher order term, 
B^ in equation (20), which is defined in analogous to B^' in equation (19) that scales 
as /jvjl- 

We consider three representative surveys: (i) idealistic survey with a huge volume 
and a deep sampling (V = 100/i~ 3 Gpc 3 , n g = 1 x 10~ 3 /i 3 Mpc~ 3 , z = 1), (ii) 
realistic survey with a large volume accessible in near future (V = 10/i _3 Gpc 3 , n g = 
5 x 10- 4 /i 3 Mpc" 3 , z = 1), and (iii) deep survey (V = 3/i" 3 Gpc 3 , n g = 3x 10- 4 /i 3 Mpc" 3 , 
z = 2). Parameters of these three surveys roughly correspond to EUCLID [37], SuMIRe 
[38] and HETDEX [39], respectively, except for the slightly smaller value of redshift 
in the deep survey. Although the mass resolution of the current simulations are not 
sufficient to reproduce the same number density of galaxies in those surveys, it is worth 
giving a rough estimate of the detectability. 

Under the assumption of the one-to-one correspondence between haloes and 
galaxies, we first estimate the minimum halo mass M min that reproduces the mean 
galaxy number density, n g , for each survey. We use the mass function of [40] to derive 
the minimum value, M min . We then compute the linear bias parameter, b\, using the 
Sheth & Tormen fit [41]. The resultant minimum masses and bias parameters are 
summarized in Tab. 1. 

Based on our numerical experiments, we focus on the isosceles triangles with 
k = 0.042/?,Mpc -1 (see the left panel of Fig. 6), and model the galaxy bispectrum 
similar to the halo bispectrum (19) as 

B g (a; z, M min ) = B<®(z, M min ) [l + f NL aC w ( Zl M min ) + / 2 L « 3 C (2) ( 2 , M min )] . (20) 

We estimate the coefficients, B^°\ and from the three halo catalogs using the 
fitting procedure described in the previous section. We then scale them to the minimum 
halo masses for the three surveys as follows: for the coefficient B^°\ we scale as cx b\. 
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This is based on the fact that the term in equation (13), which scales as bfb 2 , is the 
dominant contribution at large scales, and b 2 oc b\ at the high-peak limit (see Appendix 
B). On the other hand, we assume that and do not sensitively depend on M m i n 
and z, and treat them as constants. 

The left hand panel of Fig. 8 illustrates the scaling of measured from the N- 
body simulations. The triangles, circles and diamonds respectively correspond to the 
measurements from iV-body simulations at z — 0.5, 1 and 2, respectively. We also show 
the oc bf scaling by three solid lines. The scaling seems to be a reasonable fit to the 
simulations. In the right hand panel, we plot and as a function of the minimum 
halo mass at the three redshifts (upper. C^, lower. C^). Since we did not detect any 
significant change in these coefficients for different mass and redshift, we simply derive 
the values of and from x 2 fits to the iV-body data. We extrapolate the three 
coefficients to the halo masses corresponding to the three surveys. The accuracy of 
these scalings can only be tested using higher resolution simulations, and so leave this 
for future work. 




Figure 8. Scalings of the coefficients in equation (20). We assume that B^ oc bf, 
while and are treated as constants. We extrapolate these coefficients to the 
three surveys depicted by crosses. The triangles, circles and diamonds correspond to 
the measurements from TV-body simulations at z = 0.5, f and 2, respectively, left: B^ 
as a function of the linear bias, right: C^ 1 ' and C' 2 ' as a function of the minimum 
halo mass. 



For statistical errors of these three surveys, we consider the Gaussian contribution 
as a simple estimate, and neglect the non-Gaussian error. We have [42]: 

'P g (h) + n- 1 ] \P g {k 2 ) + n- 1 ] [P g (k 3 ) + n" 1 ] , (21) 



[AB g (h,k 2 ,k 3 )f 



-^triangle 

where iV triang i e denotes the number of independent triangular configurations in that bin, 
which roughly scales as V 2 . We count the number of triangles for each bin, A^ riangle , in 
equation (21), assuming cubic-shaped survey with the quoted volumes. See Appendix 
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10.0 



100 



Figure 9. Signal to noise ratios from the three future surveys defined in equation 
(22). Notice this is estimated from only very limited configurations of Fourier space 
triangles: isosceles with two longer sides being k\ = hi = 0.042ft,Mpc _1 . 



C, where we test this formula by comparison with TV-body simulations and show that 
it works reasonably well. We use the linear power spectrum for P g assuming the bias 
parameters listed in Tab. 1. 

Now we are in a position to discuss about the future possibility to detect the 
signature of the local-type primordial non-Gaussianity. In order to quantify the 
detectability from the a-dependence of the bispectrum, we define the signal-to-noise 
ratio: 



In evaluating Eq. (22), only the isosceles triangles with k = 0.042/iMpc -1 are used, 
and the results are shown in the right hand panel of Fig. 9. It is remarkable that 
even with the very limited number of configurations for the bispectrum, detection of 
primordial non-Gaussianity is possible in all three surveys if /nl is several dozen. This 
can be compared with the analysis neglecting the term. According to [10], using 
the full configurations of the bispectrum leads to the constraint on /nl, A/nl — 5 — 30. 
Thus, we naively expected that using the full configurations taking proper account of 
the scale dependence of bispectrum greatly enhances the detectability of primordial 
non-Gaussianity, and the constraints on /nl would be much tighter. 

Finally, it is interesting to note that the S/N for the deep survey depends steeply 
on /nl, and it exceeds that of the realistic survey of relatively large volume. This is 
primarily due to the B^> term, which scales as f^ L , being more significant at higher 
redshift, and thus helps to detect a small non-Gaussianity. In this respect, the on-going 




[B(k, a,; / NL ) - B(k, c^; / NL = 0)] 2 
AB(k,a i ;f NL = 0)* 



(22) 
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mission BOSS [43], aiming at precisely measuring the scale of baryon acoustic oscillations 
from the clustering of the LRGs at z < 0.7 and QSO absorption systems at z ~ 2.5, 
may be the promising probe for constraining or detecting primordial non-Gaussianity 
of local type. 

7. Conclusions and Discussion 

In this paper, we have studied the clustering properties of dark matter haloes 
from cosmological TV-body simulations in the presence of local-type primordial non- 
Gaussianity. We found that the halo bispectrum measured from iV-body simulations 
exhibits a strong /nl dependence which becomes most prominent for squeezed 
configurations at large scales. In particular, for realistic values of |/nl| ~ 100, the 
dependence of the halo bispectrum on /nl is well characterized by the polynomial 
expansions of /nl up to second order. Since the quadratic dependence on /nl does 
not appear in the matter bispectrum at the lowest order in perturbation theory, this 
would be a clear indicator for the existence of primordial non-Gaussianity of the local 
type. 

We have investigated the shape and scale dependence of the halo bispectrum arising 
from the f^ h term, and the simulation results are compared with theoretical predictions 
based on the local bias model. For the isosceles triangles characterized by a = ki/k 3 
and k = ki = k2, the dependence of the halo bispectrum on a measured from N- 
body simulations is found to be consistent with theoretical predictions by [23] . We also 
examined the dependence of the halo bispectrum on minimum halo mass and redshift, 
and showed that the amplitude of the halo bispectrum is more significant for more 
massive haloes at higher redshifts. 

Thus, the strong dependence of the halo/galaxy bispectrum on a makes the 
detection of primordial non-Gaussianity much more promising in future surveys. As 
a preliminary investigation, we have evaluated the signal-to-noise ratio for the scale 
dependence of the bispectrum in three representative surveys, and found that even 
with the very limited number of configurations of bispectrum it is possible to detect 
primordial non-Gaussianity if f^L is several dozen. Thus, the detectability of primordial 
non-Gaussianity is expected to be greatly improved if we use all configurations of the 
bispectrum. 

We leave the following tasks as a future work: (i) Study the effects of redshift- 
space distortions. Since we focus on very large scales, we expect that these effects 
are accurately described by linear theory, i.e., they just enhance the amplitude of the 
bispectrum in a scale independent way. (ii) Construct more elaborate theoretical models 
that are applicable to a wider range of triangular configurations and compare them with 
iV-body simulations (iii) Run higher resolution simulations where we can populate haloes 
with galaxies and measure the galaxy bispectrum directly from simulations. These tasks 
are clearly very important to exploit future surveys. 
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Appendix A. Power spectrum and bispectrum in the local bias model 

In this appendix, we compute the power spectrum and bispectrum of biased tracers 
adopting the local bias model. 

The local biasing scheme is a simple prescription to relate the galaxy/halo density 
field, <5g, to matter fluctuation, 5 m , on large scales. In this treatment, the density 
fluctuation of galaxies/haloes smoothed over the radius R, 5 g , is given by a non-linear 
function of 5 m . On large scales, it can be expanded as 

S g (x; R) = h S m (x; R) + | R) -<%} + .■ ■ (A.l) 

with or being (S^) 1 ^ 2 . For simplicity, we omit the dependence on redshifts throughout 
the appendices. Equation (A.l) can be rewritten in Fourier space as 

8 s (k;R) = b 1 5 m (k;R) 

+ \h J 7^3 [6 m (q\ R)S m (k -q;R)- (K(q; R)5 m (k - q, R))] ■ (A.2) 

Using this relation, let us consider the galaxy-matter cross spectrum. With the 
help of perturbative expansion, a straightforward calculation yields 

P gm (k; R) = hP^k; R) + h 2 J -0^B m (q, k, \k - q\;R), (A.3) 

which is valid up to third order in 5 m , and the bispectrum B m is given by the first term 
of equation (9). For the scales of our interest, the integration at the right-hand side of 
this equation can be separately done taking the large-scale limit. We then obtain [13] 

P em (t,R) =b(k,f NL ;R)P m (k;R); (A.4) 
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Here, we define 

M R (k)=M(k)W R (k), (A.6) 

with Wr being the Fourier transform of the window function. 

The above procedure can also be applied when we calculate the galaxy auto 
power spectrum, P g . However, the derivation is rather simplified if we recall that the 
deterministic bias relation holds on large scales, f Then, the auto and cross power 
spectra are tightly related with each other as {P gm (/c; R)} 2 = P g (k; R) P m (k; R), which 
leads to 

P g (k ] R) = {b(k,h L ;R)} 2 P m (k ] R) 



bo a 



2 



2 



6; t 1 + 2/NL OT)l PM (A7) 

On large scales, the effect of window functions is irrelevant, and we simply drop the 
subscript R. Introducing the bias parameter b 2 = b 2 a^, we finally obtain the galaxy 
power spectrum without smoothing: 

P s (k) = b\ jl + 2f NL ^-M~ 1 (k)^ P m (k), (A.8) 

which reproduces equation (10). In the local bias prescription, the bias parameters b± 
and b 2 are given just as the fitting parameters. On the other hand, in the halo and peak 
bias formalisms, these parameters have a specific functional form, and are related with 
each other. We will discuss this issue in Appendix B, and derive an explicit relation 
between b\ and b 2 [Eq. (B.10)]. 

Next consider the galaxy bispectrum. Again, starting from equation (A. 2), a 
straightforward calculation yields [23] 



B s (k 1 ,k 2 ,k 3 ;R) = b 3 1 



B m (h, fc 2 , k 3] R) + ^ {P m (h; R)P m (k2; R) + (eye.)} 



b 2 

+ -rB COTT (k 1 ,k 2 ,k 3 ;R) 



(A.9) 



which is valid up to fourth order in 5 m . In the above, the quantities P m and B m are the 
matter power spectrum and bispectrum, whose perturbative expressions are given in 
equations (8) and (9). On the other hand, the term -B corr represents a new contribution 
arising from the matter trispectrum, T m , and it is expressed as 

1 f d^q 

B co „{h, k 2 , k 3 ; R) = ^ j [T m {q, fci - q, k 2 , k 3 - R) + (cyc.)](A.10) 

According to the perturbative calculation by [23], the above equation can be further 
decomposed into several pieces as: 

B C oAk u k 2 , k 3 - R) = /£ L Bf {k u k 2l k 3 - R) 

J NL 

I This is indeed valid as long as we are concerned with the leading-order calculation. Sec [13] for 
alternative derivation. 



23 



+ /nl [B^ih, k 2 , k 3 ; R) + B^ih, k 2 , k 3 ; R) 



nGl, 



+AB^(k 1 ,k 2 ,k 3 ;R)J2SR(ki 



i=l 



(All) 



Here, the term .B?p comes from the leading-order trispectrum. For k £j 0.1 h Mpc it 
is approximated as 



Bf(h,k 2 ,k 3 ;R) 



2al 



x 



S-M^A^Mr^P^i) [P^(fc 2 



^(fo)] 

A<ii(?)Mii(|fei-g|)^(g) + (cyc.) 



(2tt) 3 

+ 47W i? (A; 2 )^ i? (A;3)P (A; 2 )P^(A;3 
^A4*(g)A4*(|fci - g|) 

X [P^\k 2 + q\)+P^\k 3 + q\)]- 



(eye) 



(A.12) 



with P^ being the power spectrum of <3>g- The explicit expressions for the other 
remaining terms are obtained by integrating the next-to-leading order contributions 
to the trispectrum. The resultant expressions become 



B^'(k 1 ,k 2 ,k 3 ;R) 



^W R (k l )W R (k 2 )W R (k 3 



M R (h) M R (k 2 ) 
(eye), 



B^\k u k 2 ,k 3 ;R) 



xP m (A; 1 )P m (A; 2 )P 2 w (fe 1 ,fe 2 

^WRikjWRikjP^Mik^P^h) 



(A.13) 



2al 



d 3 q y j 



WRdk, - q\)W R (q)M(\k 1 - q\)M(\k 2 + q\) 



(2tt) 3 

x [P^(\k 2 + q\) + P^ 
+ (5 permutation) 



Q 



ko , kn 



(A.14) 



The term B^® coincides with the first term of the matter bispectrum in equation (9). 
The functions J-"n(ki) and £?R,(fcj) weakly depend on the smoothing scale R. Note that 
in deriving the above expressions, we have neglected the irrelevant terms at k £ O.lh 
Mpc" 1 . 

In the expression for the galaxy bispectrum, the important findings here are a new 
term which scales as /nl and additional contributions which scale as /nl to the matter 
bispectrum. Although [23] further considered the term arising from a cubic correction, 
^nl^c to equation (1), we do not discuss about this term in this paper. See also 
[17, 22, 21] for discussions about the #nl term. Ref. [23] evaluates the asymptotic forms 
of each term at squeezed limit (a 3> 1, where k\ = k 2 = ak 3 = k) for isosceles triangular 
configurations, and the results are shown in equation (12) in the text. Note that Ref. [21] 
also compute the matter and galaxy bispectrum up to the one-loop order (i.e., 0(5q)) 
assuming local bias model for both local and equilateral type non-Gaussianity. 
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Alternatively we might be able to investigate the galaxy bispectrum using a similar 
description to the galaxy power spectrum: 

B g (fei, k 2 , k 3] R) = b(ki, / NL ; R) b{k 2 , / NL ; R) b{k 3 , / NL ; i?) S m (fci, fe 3 , fc 3 ; fl).(A.15) 

We focus on squeezed triangles at large scales, k — > and a ^> 1 and drop the window 
function. Then we get the second-order term in /nl ; which is the same as the term 
in equation (15) up to the prefactor. However, the coefficient of the term scaling as 
/nl m this prescription is 4b\ + (26/7)bfb 2 in this limit, which does not reproduce the 
contribution that depends on the smoothing scale in equation (14). This prescription 
also misses the term coming from bfb 2 [P m _(ki)P m (k2) + (cjc.)}. This implies that although 
some contributions are not included, this description captures the essence of the galaxy 
bispectrum at the squeezed limit. In other words, the /$ L term has the same origin as 
the scale dependent bias in the power spectrum. Although equation (A. 15) seems rather 
empirical, it is naturally derived in Ref. [22], where the authors calculated the galaxy 
bispectrum using a multivariate biasing scheme [see their equations (73) - (76)]. 



Appendix B. On the relation between peak biasing and local biasing models 

Scale-dependence of halo power spectrum and bispectrum in the presence of primordial 
non-Gaussianity have been derived in the literature in different ways, based on the local 
bias prescription and the halo/peak formalism. However, the resultant expressions for 
peak and halo bias coincide with each other in the high-peak/thresold limit, and there is 
a clear relationship between peak bias and local bias prescriptions. The relation between 
the linear and quadratic bias parameter plays important roles in the accurate modelling 
for the power spectrum and the bispectrum of galaxies as seen in the previous appendix. 

In this appendix, in order to elucidate these properties in a self-contained manner, 
we give an explicit relationship between the local biasing and peak biasing models, and 
show that the peak density field in the high-peak limit can be described by the local 
biasing prescription. In the end, we obtain the relation (B.10), which was used in Sec. 4 
when comparing the A-body results with model prediction of halo power spectrum. 
We also use this relation in Sec. 6 to discuss the future detectability of the local-type 
primordial non-Gaussianity. 

We begin by writing down the definition of peak density field. According to [24] , 
it is given by 

5 s (x;R) = p v (x;R) = e\6 m (x;R)-v<T R ], (B.l) 

in the Lagrangian space. In the above, is the Heaviside step function and v = 5 c /ur 
with the critical overdensity, 5 C ~ 1.686. Strictly speaking, the above definition does 
not imply the local maximum of the density field, however, the local density specified 
above is expected to roughly correspond to the peak in the high-threshold limit, v 1. 

Starting with the expression (B.l), we want to derive Taylor series expansion of 5 g 
in terms of the local density 5 m . To do this, let us first expand the peak density in terms 
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of the Hermite polynomials: 



Rn 



p v (x;R) = J2 -J H n (5 m /a K ). 



n=0 



III 



(B.2) 



The coefficient R n is given by (e.g., [44]): 

r + OO dy 



Rn 



2n 



-y 2 l 2 H n {y)Q[{y-v)a K ] 



— erfc 
2 



\V2 / 



71 = 0, 



(B.3) 



'27T 



71 > 1. 



In the high-peak limit ^ ^> 1, the coefficient i? n asymptotically approaches 

,71-1 



Rn 



-v 2 /2 



2,TT 



for n > 0. Substituting this back into (B.2), we obtain 

e -u 2 /2 oo 

p u (x;R) ~ == ^ — 7H n (6 m /a R ) 



uV2tt 



n=0 



uV2tt 



(B.4) 



(B.5) 



where we used the relation J2 n ( xn / n ^-)H n (t) = e - x2 / 2 + tx i n the last equality. Now, 
recalling from the cumulant expansion theorem, (e* t<5m ) = exp[X] n (it) n (S^ l ) c /n\], the 
averaged peak density (p u ) in the high-peak limit becomes 



(p„(x;R)) 



r,W°"R,) <5m 



vy/2ir ^\/27r 
Hence, the peak density field 5 g becomes 

v 



exp 



Ln=0 



S g (x; R) ~ exp 



Cft n 71. 



n=0 



.(B.6) 



(B.7) 



which can be expanded in the form of local biasing expression (A.l) as 



v 



1 



«5 g (a!; R) = — 5 m (x; R)+- — {5 2 Jx; R)-a 2 i }+- — {5'Jx; R)-(5 S J C }+- ■ ■ .(B.8 
cr R 2 Vcr R / 61 \a R J 

Note that this expansion is done in Lagrangian space. Assuming the usual mapping 

from Lagrangian to Eulerian space, = 1 + b^, where &e and b^ are linear bias 

parameters in Eulerian and Lagrangian space, the biasing parameters in the high-peak 

limit can be read off by comparing equation (B.8) with equation (A.l), and are expressed 

as 



u 2 u u A t z/ 6 



(B.9) 



Note that the higher-order biasing parameters should be also modified by the mapping 
from Lagrangian and Eulerian space but this effect is small in the high peaks limit and 
we simply ignore it. Then the biasing parameters have a relation 

k = 5 c (h - 1), (B.10) 

where b 2 = b 2 cr 2 l . For a better fit of halo power spectrum (10) to the N-body simulations, 
a slight modification to the above relation might be necessary [see Eq. (11)]. 
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Appendix C. The Sample Variance of the Matter and Halo Bispectrum 

It is of importance to investigate the variance of the bispectrum in the presence of 
primordial non-Gaussianity, although our simulation sets are too small to examine the 
full covariance of the bispectrum (see e.g., [45]; the authors performed 5000 realizations 
of (l/i~ 1 Gpc) 3 volume simulations to investigate it for Gaussian initial conditions). Here 
we show our measurements of the variance (i.e., the diagonal elements of the covariance 
matrix) for both the matter and halo bispectrum. 

Fig. CI shows the variance of the matter (left) and halo (right) bispectrum for the 
same configurations as in the left panel of Fig. 6 at z = 0.5. The symbols correspond to 
the variance measured from iV-body simulations, while the lines are obtained from 
equation (21). In computing equation (21), we substitute the value of the power 
spectrum measured from iV-body simulations both for matter and halo. Overall, the 
analytic predictions are good approximation of the iV-body simulations, ensuring the use 
of this formula. For the variance of the matter bispectrum, there is little evidence of /nl 
dependence. This is natural because the matter power spectrum also depends on / NL 
only weakly (see Fig. 2). On the other hand, the right panel shows a strong dependence 
on /nl, reflecting a strong /nl dependence of the halo power spectrum. The prediction 
of equation (21) seems worse at larger a and larger |/nl|- Since the leading correction 
term to this formula have the form of P(/ci)T(fc 2 , k 3 , — k 2 , — k 3 ) + (eye), this feature is 
quite reasonable. 
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Figure CI. Variance of the matter (left) and halo (right) bispectrum. The symbols 
are measured from ./V-body simulations, while the lines are computed by equation (21). 
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